## ---- Ordered logistic ----
# Note: postcodes have been pseudonomyzed
survey <- readRDS("dat/auxiliary/survey.rds")

## Estimation
m <- MASS::polr(as.factor(rent_tri) ~ rent_diff,
                data = survey,
                weights = weights,
                Hess = TRUE)
vcov_m <- sandwich::vcovCL(m, survey$postcode)

x <- seq(min(survey$rent_diff), max(survey$rent_diff), length.out = 101L)

# Parameter simulation
set.seed(20240513L)
coef_sim <- MASS::mvrnorm(10000L, c(m$coefficients, m$zeta), vcov_m)
b <- as.vector(coef_sim[, 1])
eta <- b %*% t(x)
k1 <- coef_sim[, 2]
k2 <- coef_sim[, 3]

stacked_prob <- abind::abind(t(apply(
  apply(eta, 2, function (x)
    plogis(k1 - x)), 2, quantile, c(.5, .025, .975)
)),
t(apply(
  apply(eta, 2, function (x)
    plogis(k2 - x)), 2, quantile, c(.5, .025, .975)
)),
matrix(1.0, 101L, 3L),
along = 3)

## Visualization
pdf(file = "fig/rent_perceptions.pdf",
    width = 1600 / 288,
    height = 1200 / 288)
stacked_prob_plot(
  est = stacked_prob,
  x_vals = x,
  x_lab = "Change in rents, 2015-2020",
  labels = rev(
    c("Increased\nstrongly",
      "Increased",
      "No increase/\nDecrease")
  ),
  offset = -1.5,
  rug = survey$rent_diff,
  zeroline = TRUE
)
dev.off()